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Introduction 


1.  Clinical  Relevance:  Breast  Cancer 

Breast  cancer  is  the  most  common  cancer  type  that  affects  women  globally  [11.  In  the 
United  States,  due  to  the  long  life  spans,  the  incidence  is  even  higher:  every  one  woman 
over  eight  will  develop  breast  cancer  in  her  lifetime.  It  was  estimated  that  approximately 
21 1,240  new  invasive  breast  cancer  cases  and  58,490  new  in  situ  cases  would  be  found  in 
American  women  in  2005  [2].  Moreover,  the  disease  is  the  second  leading  cause  of 
cancer-related  death  for  women  in  the  United  States,  which  was  predicted  to  kill  40,870 
women  in  2005  [2], 

Presently  there  is  no  effective  way  of  preventing  the  disease.  However,  the  detection  of 
the  cancer  at  its  early  stage  has  been  found  to  significantly  improve  the  survival  rates  [3- 
61.  For  example,  when  the  breast  cancer  is  detected  at  the  localized  stage,  the  five  year 
relative  survival  rate  is  98%  [2],  By  contrast,  when  it  is  not  found  until  metastasized,  the 
five  year  survival  rate  drops  dramatically.  In  addition,  when  the  cancer  is  found  earlier, 
more  viable  treatment  options  are  also  available  [7-91. 

X-ray  mammography  is  a  successful  tool  for  the  early  detection  of  the  breast  cancer.  The 
abnormalities  can  manifest  themselves  on  a  mammogram  as  either  masses,  clusters  of 
microcalcifications,  or  architectural  distortions  even  before  any  symptom  shows  up.  An 
annual  screening  program  based  on  mammography  is  recommended  for  women  older 
than  forty  years  or  younger  women  with  high  risk  by  National  Cancer  Institute,  American 
Cancer  Society  and  American  College  of  Radiology.  It  has  been  proven  to  reduce  the 
mortality  rate  of  the  breast  cancer  since  its  initiation.  For  example,  screening 
mammography  decreases  the  fifteen  year  mortality  for  women  in  their  forties  by  20% 
[101.  Also,  it  is  found  that  screening  is  most  effective  for  women  older  than  55  years  old 
[10,  111. 

2.  Limitation  of  Screening  Mammography 

Film-screen  X-ray  mammography  is  presently  the  only  FDA  approved  screening  tool 
aiming  at  early  detection  of  the  breast  cancer.  While  it  has  been  proven  to  be  effective,  it 
is  not  omnipotent  in  its  detection  sensitivity  of  breast  lesions  due  to  several  limitations 
such  as  two-dimensional  (2D)  projection  data  acquisition  and  restricted  range  of  linear 
optical  response  of  the  detector.  Overall,  it  has  a  sensitivity  within  the  range  of  63%  to 
88%  depending  on  the  patient’s  age  group,  family  history  [121  and  breast  density  [131. 
For  women  with  dense  breasts,  the  sensitivity  is  lower  since  in  their  mammograms  the 
dense  appearance  of  the  breast  tissue  is  more  likely  to  obscure  any  abnormalities  and 
makes  the  detection  of  breast  cancer  even  more  challenging  [141.  In  addition,  the 
situation  gets  complicated  by  the  fact  that  breast  density  is  also  a  risk  factor,  which  means 
that  women  with  dense  breasts  tend  to  be  more  likely  to  get  breast  cancer. 

3 .  Emerging  Dedicated  Breast  CT  Imaging 

Breast  CT  technology  offers  the  potential  to  detect  breast  lesions  among  women  with 
dense  breasts,  which  is  a  well-known  challenging  task  in  conventional  film-screen  X-ray 
mammography  [151. 
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The  breast  imaging  via  CT  modality  was  first  conducted  in  the  late  1970s  [16,17].  At  that 
time  the  breast  images  were  acquired  slice-by-slice  using  x-ray  fan-beam  geometry.  The 
long  scan  time,  high  patient  dose  together  with  the  limited  image  quality  prevented  the 
further  application  of  this  technique. 

With  the  development  of  flat-panel  detectors  in  recent  years,  the  fast  breast  imaging  via  a 
cone-beam  CT  became  possible  and  breast  CT  regained  attention.  So  far,  four  research 
groups  have  fabricated  their  own  dedicated  breast  CT  systems  and  one  research  group  has 
set  up  a  virtual  breast  CT  system  by  applying  mathematical  models.  The  four  groups  with 
actual  breast  CT  systems  are:  Dr.  Boone’s  group  in  University  of  California,  Davis 
[18,19],  Dr.  Tornai’s  group  in  Duke  University  [20,21],  Dr.  Ning’s  group  in  University  of 
Rochester  [22,23]  and  Dr.  Shaw’s  group  in  University  of  Texas  M.D.  Anderson  Cancer 
Center  [24,  25].  The  one  group  with  a  virtual  breast  CT  system  is  Dr.  Glick’s  group  in 
University  of  Massachusetts  [26,  27]. 

In  a  conventional  CT  system,  the  x-ray  tube/  detector  move  around  the  torso  of  a  patient; 
whereas  a  dedicated  breast  CT  system  has  a  joint  x-ray  tube/detector  move  just  around 
the  breast.  The  system  is  normally  set  up  in  the  following  way:  a  patient  lies  supinely  on 
a  table  with  one  breast  hanging  through  a  hole.  The  flat-panel  detector  is  installed 
vertically  and  moves  jointly  with  the  x-ray  tube. 

By  this  design  of  the  dedicated  system,  the  field  of  view  (FOV)  of  the  detector  can  be 
fully  employed  for  breast  imaging.  What’s  more,  since  no  other  tissues  will  attenuate  the 
x-ray  beam,  the  effective  glandular  dose  delivered  to  the  patient  can  be  lowered  to  match 
the  two-view  screening  mammograms  for  the  same  breast  [28]. 

Although  they  share  the  same  physical  setup,  the  four  breast  CT  systems  differ  in  their 
detailed  technical  aspects:  the  choice  of  x-ray  beam,  the  x-ray  source  orbit,  and  the  peak 
voltage  and  tube  current  values  used. 

The  proposed  project  is  a  collaborative  effort  between  our  group  at  Duke  and  Dr. 
Boone’s  group  in  University  of  California,  Davis.  Based  on  the  raw  data  provided  by 
them,  we  will  develop  the  techniques  for  image  improvement  via  the  scatter 
compensation  and/or  denoising. 


Report  Body 

In  the  approved  statement  of  work  (SOW),  it  was  proposed  that  task  1  to  3  would  be 
finished  within  the  first  fiscal  year.  What  has  been  done  so  far  basically  follows  the  SOW 
with  some  minor  changes,  which  will  be  stated  in  each  subsection. 

Task  1:  Develop  and  test  a  unique  two-dimensional  Bayesian  image  processing  technique 
on  the  projection  data  of  cone-beam  breast  Computed  Tomography  ( breast  CT)  obtained 
without  a  grid. 

In  the  previous  work  [29],  the  Bayesian  image  estimation  technique  based  on  Poisson 
noise  model  has  been  successfully  applied  to  reduce  the  scattering  radiation  in  digital 
mammograms.  While  CT  projections  of  the  breast  are  similar  enough  to  give  us 
confidence  in  its  eventual  success,  this  application  of  BIP  is  unique  and  will  present 
novel  challenges.  For  example,  in  digital  mammography,  the  projection  images  are  the 
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ultimate  images  for  display  and  diagnosis;  whereas  in  breast  CT,  the  diagnosis  is  based 
on  the  reconstructed  images  or  its  derivatives  such  as  maximum  intensity  projection  (MIP) 
[301.  To  reconstruct  the  tomographic  images,  the  projection  images  need  to  be  converted 
the  line  integrals  of  attenuation  coefficient  via  the  logarithmic  operation.  It  is  expected 
the  characteristics  of  the  noise  will  be  different  because  of  the  nonlinear  operation. 
Therefore  the  component  for  noise  constraining  in  the  BIP,  i.e.  Gibbs  prior,  has  been 
modified  to  incorporate  this  knowledge. 

In  addition,  while  in  previous  implementations  of  BIP  for  digital  chest  and  breast  imaging 
a  Poisson  noise  model  was  used,  the  flat-panel  detectors  actually  are  integrating  detectors 
and  so  the  Poisson  noise  model  is  not  exactly  correct.  A  more  exact  statistical  model  for 
the  energy-integrating  detector  has  been  incorporated  into  the  scatter  compensation 
framework  and  the  maximum  likelihood  estimate  (MLE)  of  scatter-free  image  has  been 
acquired  via  the  expectation  maximization  (EM)  algorithm.  The  details  are  shown  as 
follows. 


1)  Signal  from  an  Energy-Integrating  Detector 

To  simplify  the  notion,  assume  that  the  x-ray  spectrum  impinging  on  the  energy- 
integrating  detector  is  binned  into  the  discrete  energy  bands.  In  our  analysis,  the  primary 
x-ray  photons  surviving  attenuation  and  the  scattered  photons  will  be  analyzed  separately, 
so  there  are  x-ray  spectrums  for  primary  and  scatter  radiations  respectively.  Moreover, 
the  x-ray  spectrums  are  dependent  on  the  attenuation  material  and  its  thickness. 


n  ~  Poisson(A) 

e  =  K±ei 

i= 1 


(1) 


where  n  is  the  number  of  photons  hitting  on  the  detector,  ej  is  the  energy  carried  by 
individual  photon  i  and  follows  probability  outlined  by  the  x-ray  spectrum,  K  is  the  gain 
factor  and  e  is  the  signal  recorded  by  the  energy-integrating  detector. 

Obviously,  the  signal  e  does  not  follow  a  Poisson  distribution,  but  rather  a  compound 
Poisson  distribution  [311.  Assuming  the  first  and  second  central  moments  of  x-ray 
spectrum  for  e,  are  p  and  a",  because  of  the  mutual  independence  between  ej’s,  we  get  the 
first  and  second  central  moments  of  e  as: 

E(e)  =  e  =  En(E(Kfje,\  «))  =  En  (K^  E(e,  I  n)  =  En  (. Knju )  =  KA/u  ,  (2) 

i= 1  i= 1 

n  n 

var(e)  =  varn(£'(A'^e;  I n))  +  En(var(K'^jej  I n))  =  var n(Kn/j)  +  En(K2ncr2) 

1=1  1=1 

=  K2(Afi2  +  A<j2)  =  e  ■  K—  +a  •  (3) 

P 

Therefore,  a  Gaussian  approximation  of  e  is: 
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e  ~  Gaussian{e,e  K—  +<T  ) 


where  ( )+  means  that  e  can  only  take  on  the  non-negative  values. 

2)  Gaussian  noise  model 

The  Gaussian  noise  model  for  the  scatter  radiation  is  proposed  as  follows: 

2  2 
di  I  B. cr(1  ~  Gaussian(bi,<JiX  ) 

y  I  B,crn2  ~  Gaussian((B**P)j,cri22) 

yt  =  d)  +  Sj  I  B,crn2 ,<jj22  -  Gaussian (b,  +  ( B**P)i,crn "  +  cr22) 


where  di,  s;,  and  y,  represent  respectively  the  primary,  scatter  and  total  radiations  reaching 
a  pixel  i.  The  b;  is  the  expectation  of  the  primary  radiation  di  at  pixel  i.  B  ={  b; ,  i=l,...,N} 
is  its  two-dimensional  matrix  representation.  In  addition,  o,f  and  Oj2‘  represent  the 
variance  of  the  primary  radiation  and  the  variance  of  the  scatter  radiation  in  each  pixel  i. 

The  EM  algorithm  for  the  MLE  estimate  of  B  based  on  the  Gaussian  noise  model  is 
derived  via  the  steps  shown  in  Appendix  2,  and  the  updating  equation  for  bk  at  iteration 
step  n+1  is: 

bt"  =  +wk-[yk- (V'0  +  (5W  * */*)*)] ;  (5) 

where 

wk~(7ki  bcrkl  +<yk2  )  (61 


Incorporating  the  equation  (3)  into  here,  we  estimate  Wk  in  an  iterative  manner: 


w-bl  w-Mp  +C7p  /(A (>,)  •  — 


2  2  2  2 
^ P  |  PjL  ^ s  ^ 

MP  "  Ms 


=  bkn)  /(bkw  +  ( B(n)  *  *P),  ■  (A  +a'  /Mp  +(7p  )) 


Other  than  the  existence  of  a  nice  analytic  formula  for  b,  approximating  the  sum  of 
energies  e  by  a  Gaussian  distribution  poses  an  additional  advantage:  the  dark  current  (dc) 
and  electronic  readout  noise  can  be  easily  incorporated  into  the  framework.  The 
electronic  readout  noise  is  usually  represented  by  a  Gaussian  distribution  with  variance  of 
Gem”-  Then  both  components  are  jointly  modeled  by  N(dc,oern  )  [32].  Incorporating  this 
into  the  model  shown  in  equation  (4),  the  signal  is  now: 

2  _j_  2 

e  ~  Gaussian{e  +  dc,e  ■  K  — - —  +  cr  2)+ .  (8) 
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Task  2:  Reconstruct  the  three-dimensional  breast  image  based  on  the  processed 
projection  data  from  Task  1. 

An  existing  ordered  subset  expectation  maximization  (OSEM)  algorithm  for  cone-beam 
CT  was  configured  for  tomographic  reconstruction.  However,  there  was  some  severe 
artifact  in  the  reconstructed  volume.  The  cause  of  the  artifact  has  not  been  pinpointed  yet. 
As  an  alternative,  the  Feldkamp  type  filtered  back-projection  (FBP)  algorithm  [33]  was 
implemented.  Fortunately,  due  to  the  sufficient  sampling  of  the  projection  data  and 
projection  views,  the  FBP  algorithm  provided  good  reconstruction  results.  It  is  also 
computationally  more  efficient  than  the  OSEM  algorithm.  Therefore,  we  feel  that  FBP 
algorithm  is  a  good  substitute  for  the  OSEM  algorithm  for  Task  2  in  the  approved  SOW. 

Task  3:  Apply  the  algorithm  in  Task  1  to  the  two-dimensional  slices  of  the  reconstructed 
three-dimensional  breast  image  from  the  unprocessed  projection  data. 

The  BIP  algorithm  in  Task  1  has  been  set  up  primarily  for  image  noise  reduction  and  tried 
on  the  reconstructed  slice  image.  Visually,  the  processed  images  are  smoother  than  the 
original  ones.  Further  evaluation  will  be  conducted. 

Within  the  next  two  fiscal  years,  task  4  to  6  will  be  fulfilled. 

Task  4:  Develop  and  test  three-dimensional  Bayesian  image  processing  technique  on  the 
reconstructed  image  based  on  the  unprocessed  projection  data  acquired  withou  t  a  grid. 
Task  5:  Develop  a  Computer  Aided  Diagnosis  tool  for  detecting  breast  mass  lesions 
based  on  the  projection  data. 

Task  6:  Test  and  compare  the  performances  of  the  CAD  developed  in  Task  5  applied  to 
processed  projection  data  from  Task  1  with  the  CAD  performance  on  the  projection  data 
without  Bayesian  processing. 

Due  to  the  limited  cases  available,  a  major  change  to  Task  5  will  be  an  addition  of  3D 
mass  simulation  routine  to  generate  multiple  synthetic  datasets  as  well  as  simplified  ROC 
analysis.  It  is  felt  that  this  change  won’t  affect  the  overall  structure  of  the  project,  since 
the  primary  goal  of  the  proposed  CAD  tool  is  for  the  evaluation  of  the  image  processing 
techniques  developed  during  the  first  two  years  of  the  grant. 

Key  Research  Accomplishments 

•  Proposed  a  Gaussian  noise  model  for  scatter  compensation  when  an  energy- 
integrating  detector  is  used  for  image  acquisition; 

•  Derived  the  EM  algorithm  for  the  MFE  estimate  based  on  the  Gaussian  noise 
model; 

•  Implemented  the  Feldkamp  type  FBP  for  the  cone-beam  CT; 

•  Reconstructed  the  breast  CT  data; 

•  Implemented  a  partial  diffusion  equation  (PDE)  based  image  denoising  technique; 

•  Compare  the  reconstructed  volumes  with  and  without  the  image  processing  by  the 
metrics  such  as  CNR; 

•  Coded  the  Siddon  algorithm  for  forward  projection  of  3D  voxelized  data  using 
cone-beam  geometry. 
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Reportable  Outcomes 


•  The  PI  finished  a  project  titled  ‘Gaussian  Noise  Model  for  Scatter  Compensation 
in  Digital  Mammography’  and  earned  her  Master  of  Science  degree  in  statistics 
from  Institute  of  Statistics  and  Decision  Science  (ISDS),  Duke  University.  Please 
refer  to  Appendix  2. 

•  The  PI  applied  for  the  Sally  Hughes-Schrader  travel  grant  of  Duke  University  for 
a  potential  site  visit  to  the  breast  CT  research  group  in  University  of  California, 
Davis  in  the  summer  of  2006. 

Conclusions 

In  summary,  the  development  of  dedicated  cone-beam  breast  CT  imaging  technique 
requires  effective  and  efficient  ways  to  reduce  its  scattered  radiation  as  well  as  denoising 
so  as  to  provide  high  quality  images  and  help  make  diagnostic  decisions.  Therefore,  it  is 
important  to  investigate  different  possible  image  processing  tools  and  decide  which  one  is 
better  based  on  image  quality  metrics  such  as  CNR  as  well  as  the  observer  performance 
study  via  a  ROC  analysis. 

A  Gaussian  noise  model  accounting  for  the  energy-integrating  characteristic  of  a  flat- 
panel  detector  has  been  developed.  Based  on  this  model,  the  MLE  estimate  of  the  scatter 
free  image  via  the  EM  algorithm  has  been  derived.  A  partial  diffusion  equation  based 
image  denoising  technique  has  been  implemented. 

The  raw  breast  CT  data  has  been  successfully  reconstructed  via  the  Feldkamp  FBP 
algorithm  for  cone-beam  geometry.  The  reconstructed  images  totally  eliminate  the 
problem  caused  by  the  superposition  of  breast  tissues  which  is  the  limitation  of  the 
conventional  screening  mammography. 
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Abstract 


Breast  cancer  is  the  most  common  cancer  type  that  affects  women  worldwide.  In  the 
United  States,  every  one  woman  over  eight  will  develop  breast  cancer  in  her  lifetime. 
Although  no  effective  way  of  preventing  the  disease  has  been  found,  early  detection 
of  the  cancer  through  noninvasive  breast  imaging  is  desirable  because  it  warrants 
more  choices  of  viable  treatments  and  higher  survival  rates.  Digital  Mammography  is 
among  such  imaging  techniques. 

Compton  scattering  of  x-ray  photons  is  one  mechanism  of  attenuating  the  x-ray 
beam,  which  in  turn  forms  the  contrast  in  a  projection  image.  However,  its  detection 
in  the  projection  image  is  a  cause  of  image  quality  degradation  since  it  will  add  noise 
to  the  image  and  reduce  the  contrast.  Therefore  many  efforts  are  made  to  reduce  the 
detected  scatter  radiation  in  the  projection  image  either  by  applying  some  hardware 
during  acquisition  or  by  using  post-acquisition  software  compensation.  The  method 
presented  in  this  thesis  belongs  to  the  latter  category.  A  Gaussian  noise  model  for 
scatter  is  proposed  and  its  EM  ML  estimation  is  derived.  In  addition,  Bayesian  MAP 
estimation  is  obtained  by  applying  a  Gibbs  prior  with  a  discontinuity  adaptive 
potential  function. 

The  previously  proposed  Poisson  noise  model  is  flawed  in  that  the  radiation  does 
not  directly  follow  a  Poisson  distribution.  Instead,  a  Gaussian  distribution  can 
reasonably  describe  the  radiation  data.  When  a  computation  method  like  Gibbs 


sampling  is  used,  Poisson  noise  model  will  give  erroneous  results  due  to  the  incorrect 
modeling.  The  conversion  factor  between  radiation  and  the  number  of  photons  is 
energy  dependent.  If  it  is  approximated  by  a  constant  independent  of  energy,  then  the 
Poisson  noise  model  can  be  justified  through  a  latent  data  augmentation  scheme 
when  EM  algorithm  is  used. 

The  reason  we  pursue  the  EM  computation  in  this  thesis  is  that  it  has  a  nice 
analytic  formula.  Due  to  the  large  number  of  pixels  in  an  image  and  the  existence  of 
a  convolution  operation,  the  computation  can  be  greatly  reduced  with  this  analytic 
formula. 

The  digital  mammography  image  of  a  uniform  breast  phantom  is  processed  by 
the  MLE  and  MAP  algorithms  of  Poisson  model  and  Gaussian  model.  The  results  are 
compared  through  the  image  quality  metrics  like  the  residual  scatter  radiation,  the 
contrast- to-noise  ratio  and  the  spatial  resolution. 

From  the  results  we  get,  it  is  shown  that  Gaussian  noise  model  can  be  used  to 
reduce  the  scatter  radiation  in  the  digital  mammography  images.  Its  performance  is 
improved  by  incorporating  Gibbs  priors  without  loss  of  resolution.  In  addition, 
Gaussian  noise  model  works  slightly  better  in  improving  the  contrast-to-noise  ratio 


than  the  Poisson  noise  model. 
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Chapter  1 
Introduction 


1 . 1  Digital  Mammography 

Breast  cancer  is  the  most  common  cancer  type  that  affects  women  worldwide  [  1  ] .  In 
the  United  States,  every  one  woman  over  eight  will  develop  breast  cancer  in  her 
lifetime.  And  it  was  estimated  that  approximately  211,240  new  cases  of  invasive 
breast  cancer  would  be  found  in  American  women  and  the  disease  would  kill  40,410 
women  in  2005  [2] .  Although  there  is  no  effective  way  of  preventing  the  disease,  it  is 
desirable  to  find  early  signs  of  the  cancer  (e.g.,  impalpable  masses  and/or 
micro-calcifications)  through  noninvasive  breast  imaging  techniques  such  as  x-ray 
mammography.  The  detection  of  the  cancer  at  its  early  stage  warrants  more  choices 
of  viable  treatments  and  higher  survival  rates[3-51. 

An  x-ray  mammography  system  is  typically  comprised  of  two  major  parts:  the 
x-ray  source  and  the  detector.  Depending  on  the  type  of  the  detector,  the 
mammography  system  can  be  either  analog  or  digital.  The  analog  system  utilizes  a 
screen-film  as  the  detector,  and  is  the  only  FDA  approved  screening  tool  aiming  at 
the  early  detection  of  the  breast  cancer  for  women  more  than  40  years  old.  While  it 
has  been  proven  to  be  effective,  it  has  several  shortcomings:  1)  the  analog  film  has 
narrow  latitude.  Overexposure  or  underexposure  of  the  film  will  result  in  a  poor 
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image  which  will  be  unacceptable  for  breast  cancer  detection  and  diagnosis.  2)  The 
film  development  is  critical  for  the  quality  of  the  mammogram  as  well  as 
time-consuming.  3)  The  radiology  department  needs  a  lot  of  space  and  personnel  to 
keep  the  films.  And  4)  the  transfer  of  the  films  between  departments  or  hospitals  not 
only  is  a  lot  of  hassle,  but  also  causes  the  wear  and  tear  of  films,  which  is  inevitable 
since  they  are  often  the  only  copies  of  the  case.  Due  to  these  limitations  of  using 
films  as  the  recording  media,  many  radiology  departments  are  trying  to  go  film-less. 
It  is  realized  by  the  digital  mammography  technique. 

A  digital  mammography  system  utilizes  a  flat-panel  detector  instead  of  a 
screen-film  detector.  Recent  studies  show  that  the  diagnostic  accuracy  based  on 
digital  mammograms  is  comparable  to  those  based  on  conventional  film 
mammograms  [6,  7].  In  some  situations,  the  digital  mammography  works  even  better 
[6].  In  addition,  a  digital  mammography  system  enjoys  the  following  merits.  After 
x-ray  exposure,  a  digital  image  of  the  breast  can  be  readily  read  out  from  the 
flat-panel  detector  within  seconds.  There  is  no  overexposure  or  underexposure  issue 
related  with  this  type  of  image  since  the  flat-panel  detectors  have  a  wide  latitude  and 
excellent  linear  relationship  between  pixel  values  and  exposure  levels.  The  image 
can  be  saved  into  different  media.  The  transfer  and  copy  of  the  images  are  easy,  fast 
and  reliable.  A  very  recent  study  [81  shows  that  the  digital  mammography  images  can 
be  accurately  transferred  via  the  broadband  internet,  which  will  greatly  alleviate  the 
problem  related  to  the  shortage  of  mammographers  as  well  as  improve  the  accuracy 
of  the  diagnosis.  Moreover,  the  digital  format  of  images  makes  advanced  imaging 
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(e.g.,  breast  tomosynthesis  [9,  10])  and  image  processing  techniques  (such  as  the 


technique  presented  in  the  thesis)  feasible. 


1.2  Scatter  Radiation  and  Its  Degrading  Effect  on  the  Quality  of  Medical 
Images 


X-ray  source  emits  x-ray  photons  with  different  energies.  Figure  1.1  shows  a 


Figure  1.1:  Typical  energy  spectrum  of  an  x-ray  beam.  The  abscissa  represents  the 
energy  levels  that  a  photon  can  possibly  take  on,  and  ordinate  represents  the  number 
of  photons  having  the  corresponding  energy  level.  The  peak  voltage  (in  the  unit  of 
kVp)  corresponds  to  maximal  energy  level.  The  spectrum  shown  here  has  a  peak 
voltage  of  52.5  kVp. 


typical  energy  spectrum  of  an  x-ray  beam.  When  the  beam  passes  through  the  object 
to  be  imaged  such  as  a  breast,  it  interacts  with  the  matter  and  gets  attenuated.  At  the 
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diagnostic  energy  level,  there  are  three  basic  mechanisms  for  x-ray  attenuation: 
photoelectric  effect,  Compton  scattering  and  Raleigh  scattering  [11].  Raleigh 
scattering,  also  called  coherent  scattering,  accounts  for  less  than  5%  of  the  total 
interactions  between  x-rays  and  the  matter.  Therefore  it  is  often  omitted  for 
consideration.  Photoelectric  effect  occurs  when  x-ray  photons  are  totally  absorbed  by 
the  atoms  within  the  tissue,  as  illustrated  by  ray  1  in  Figure  1.2.  Compton  scattering 
occurs  when  the  photons  are  deflected  from  their  incident  path  with  partial  energy 
loss.  These  photons  are  called  scattered  photons  or  scatter  radiation,  as  shown  by  ray 
2  in  Figure  1.2.  The  rest  will  survive  the  attenuation  and  are  called  the  primary 
photons  or  primary  radiation  (ray  3  in  Figure  1.2).  The  primary  radiation  differs  from 
location  to  location,  which  forms  the  contrast  of  various  tissues  in  the  final  image. 

The  scatter  radiation  escaped  from  the  imaged  object  can  either  miss  the  detector 
or  impinge  on  it.  The  latter  will  be  inevitably  detected  due  to  the  fact  that  the  detector 
typically  has  broad  energy  sensitivity  and  does  not  effectively  reject  photons  that 
have  lost  energy  by  scattering.  The  total  radiation  detected  is  thus  the  sum  of  the 
primary  radiation  and  the  scatter  radiation.  The  detection  of  scattered  photons  in 
locations  that  are  different  from  their  original  path  will  add  a  component  of  noise  and 
cause  the  blurring  of  the  image. 
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X-ray  Beam 


Figure  1.2:  Illustration  of  possible  interactions  between  x-ray  photons  and  the  matter 
within  the  diagnostic  x-ray  energy  range.  The  ellipsoid  represents  the  object  to  be 
imaged.  Photons  can  be  totally  absorbed  by  the  photoelectric  effect  (ray  1),  or  be 
scattered  through  Compton  scattering  (ray  2)  and  Raleigh  scattering  (a  very  small 
portion,  thus  neglected).  The  rest  will  survive  the  attenuation  and  are  called  the 
primary  photons  or  primary  radiation  (ray  3). 


To  see  how  detected  scatter  radiation  adversely  affect  the  quality  of  a  medical 
image,  let’s  take  a  look  at  the  simple  example  shown  in  Figure  1.3.  An  ellipsoid 
lesion  is  embedded  in  a  uniform  background.  In  the  ideal  case  where  no  scatter 
radiation  is  detected,  the  total  radiation  is  equal  to  the  primary  radiation.  In  Figure 
1.3  (a),  assume  radiation  in  the  background  is  Po  and  radiation  in  the  lesion  is  Pi, 
then  the  contrast  of  the  lesion  is  Ci  =  (Pi-Po)/Po.  If  in  practice,  a  constant  scatter 
radiation  S  is  added  all  over  the  image  (shown  in  Figure  1.3(b)),  then  the  radiation  in 
the  background  becomes  To=  Po+S  and  radiation  in  the  lesion  becomes  T  i  =  Pi+S. 
The  contrast  of  the  lesion  in  this  image  will  be  C2  =  (Ti-To)/To=  (P1-P0)  /  (Po+S)  = 
Ci  *[P0/  (Po+S)],  which  is  smaller  than  Ci.  That  is,  the  detection  of  scatter  radiation 
reduces  the  contrast  of  the  lesion.  In  cases  where  scatter  is  large,  lesions  can  even  be 
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obscured. 


(a)  ideal  image  with-aul  scalier 


Figure  1.3:  A  simple  example  to  demonstrate  the  adverse  effect  of  detected  scatter 
radiation  on  image  quality.  An  ellipsoid  lesion  is  embedded  in  a  uniform  background, 
(a)  The  ideal  image  without  scatter  radiation  is  shown,  (b)  The  actual  image  is  shown, 
which  is  the  sum  of  primary  radiation  and  scatter  radiation.  The  contrast  of  the  lesion 
is  decreased. 


1.3  Scatter  Compensation 

In  summary,  the  scatter  radiation  is  a  physical  phenomenon,  which  together  with 
photoelectric  effect  causes  the  attenuation  of  the  x-ray  beam.  The  detection  of  scatter 
radiation  on  the  detector  will  degrade  the  quality  of  the  image  and  thus  adversely 
affect  the  medical  diagnosis.  This  issue  exists  widely  in  many  imaging  techniques 
such  as  Single  Photon  Emission  Computed  Tomography  (SPECT)  [12],  Positron 
Emission  Tomography  (PET)  [13],  and  projection  radiographies  like  chest 
radiography  [14]  and  mammography  [151.  Therefore,  it  is  important  to  reduce  the 
scatter  radiation  that  is  detected.  Or  in  other  words,  scatter  radiation  needs  to  be 
compensated. 

There  are  two  general  categories  of  scatter  radiation  compensation  methods:  one 
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is  hardware  compensation  such  as  the  application  of  anti-scatter  grids  [161,  slot 
scanning  systems  [17],  or  air  gaps  [18];  the  other  is  software  compensation  via 
post-acquisition  image  processing,  such  as  simple  estimation-subtraction  [19], 
convolution-subtraction  [20],  de-convolution  [21],  artificial  neural  networks  [22], 
maximum  likelihood  expectation  maximization  (EM-MLE)  [23],  or  Bayesian  image 
estimation  [24,  251. 


X-ray  Beam 


Figure  1.4:  An  anti-scatter  grid  can  be  added  on  top  of  the  detector  to  remove  scatter 
radiation.  However,  as  shown  by  the  middle  ray,  some  primary  radiation  will  be 
blocked  as  well.  To  maintain  the  image  quality,  the  patient  dose  has  to  be  increased. 


An  anti- scatter  grid  is  routinely  used  on  a  clinical  screen-film  mammography 
system.  Figure  1.4  illustrates  how  an  anti-scatter  grid  can  be  used  to  reduce  the 
scatter  radiation.  The  orientation  of  the  grid  slots  is  parallel  to  the  primary  radiation. 
Most  primary  radiation  will  pass  through  the  slots  and  reach  the  detector.  Scatter 
radiation,  by  contrast,  will  mostly  hit  on  the  metal  slits  and  be  absorbed  by  them. 
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Thus  an  anti-scatter  grid  effectively  removes  many  scatter  radiation.  Its  major 
drawback  is  that  it  also  removes  some  primary  radiation,  as  shown  by  the  middle  ray 
of  Figure  1.4.  To  maintain  the  same  image  quality,  the  magnitude  of  the  x-ray  beam 
needs  to  be  increased,  which  will  also  increase  the  total  absorbed  dose  of  the  patient. 

By  contrast,  post-acquisition  image  processing  techniques  won’t  change  the  dose 
that  a  patient  receives.  In  addition,  some  studies  [15,  24]  show  that  they  can  be  more 
effective  than  an  anti-scatter  grid  in  scatter  compensation. 

A  fundamental  assumption  behind  the  image  processing  techniques  is  that  the 
scatter  radiation  can  be  approximated  by  the  convolution  of  the  primary  radiation  and 
a  scatter  kernel.  It  is  verified  both  theoretically  [261  and  empirically  [27,  28].  In  the 
two-dimensional  case,  if  Y  is  used  to  represent  the  matrix  of  detected  total  radiation 
at  each  pixel,  D  for  the  matrix  of  the  primary  radiation,  S  for  the  matrix  of  the  scatter 
radiation,  and  P  for  the  matrix  of  the  scatter  kernel,  then  the  following  equation  is 
true: 

Y  =  D  +  S  =  D  +  D**P  =  D**(S  +  P),  (1.1) 

where  **  is  the  two-dimensional  convolution  operator  and  5  is  the  Dirac  delta 
function  in  a  matrix  form.  The  task  of  scatter  compensation  is  equivalent  to 
estimating  the  unknown  D  from  the  measured  Y. 

One  solution  is  to  de-convolve  the  Equation  (1.1)  [21,  29,  30].  If  this  is  done 
through  the  Fourier  Transform  (FT),  then 

D  =  FT  '(  FT(Y)  ).  (1.2) 

FT{S+P) 
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Or,  statistical  models  can  be  formulated  to  solve  the  problem.  In  the  past  Poisson 
noise  model  was  used,  which  assumed  that  both  primary  radiation  and  scatter 
radiation  follow  Poisson  distributions.  The  maximum  likelihood  estimate  (MLE)  of 
D  in  two-dimensional  projection  radiography  was  obtained  by  borrowing  the 
iterative  equation  originally  derived  for  SPECT  reconstruction  [23],  and  it  was 
combined  with  a  Gibbs  prior  to  form  the  maximum  a  posteriori  (MAP)  estimator  of 
D  [24].  Later  a  revision  was  made  on  the  iterative  equation  [26].  Although  promising, 
the  Poisson  noise  model  presents  some  problems:  1)  the  primary  radiation  and  scatter 
radiation  can  not  be  directly  modeled  to  follow  Poisson  distribution;  2)  due  to  the 
polychromatic  characteristic  of  x-ray  beam,  the  radiation  is  not  only  related  to  the 
number  of  photons  but  to  the  energies  of  the  photons  as  well.  Therefore,  in  this  thesis, 
a  new  explanation  is  given  that  justifies  the  Poisson  noise  model,  and  a  new  model  is 
proposed,  implemented  and  tested  on  the  digital  mammography  data  for  the 
reduction  of  detected  scatter  radiation. 

1 .4  Overview  of  the  Thesis 

The  thesis  is  organized  as  follows.  In  Chapter  2,  the  old  Poisson  noise  model  is 
briefly  introduced.  The  Gaussian  noise  model  is  then  proposed  and  its  analytical  EM 
algorithm  is  derived.  Moreover,  the  Gibbs  prior  is  incorporated  into  the  algorithm  to 
constrain  the  noise  in  the  processed  image.  Chapter  3  presents  the  latent  data 
augmentation  scheme  to  justify  the  Poisson  noise  model,  the  image  processing 
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results  obtained  from  both  the  Poisson  model  and  the  Gaussian  model,  their 


comparison  as  well  as  the  some  further  evaluation  of  the  Gaussian  noise  model.  The 
thesis  is  concluded  in  chapter  5. 
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Chapter  2 

Materials  and  Methods 


2.1  Scatter  Kernel 

Towards  the  end  of  the  previous  chapter,  we  mentioned  that  scatter  radiation  is  often 
modeled  as  the  convolution  of  the  primary  radiation  and  a  scatter  kernel.  The  scatter 
kernel  is  also  called  scatter  point  spread  function  (PSF).  Experiments  [31,  32]  and 
Monte-Carlo  simulation  models  [33,  34]  showed  that  the  scatter  kernel  can  be 
represented  by  a  circularly  symmetric  exponential  decay  curve.  The  curve  can  be 
uniquely  determined  by  two  parameters:  the  magnitude  (M)  and  full  width  at  half 
maximum  (FWHM).  The  two  parameters  are  illustrated  in  Figure  2.1  (a).  Figure 
2.1(b)  shows  the  three-dimensional  representation  of  such  a  scatter  kernel  with  M  of 
1  and  FWHM  of  80  pixels.  Throughout  the  thesis,  matrix  P  is  used  to  represent  the 
scatter  kernel  and  it  is  known  a  priori. 
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Figure  2.1:  A  sample  scatter  kernel  with  magnitude  of  1  and  FWHM  of  80  pixels,  (a) 
One-dimensional  profile  of  the  kernel,  (b)  Its  three-dimensional  surface  plot. 


2.2  Poisson  Noise  Model 

When  put  in  a  statistical  framework,  equation  (1.1)  becomes: 

E(Y)  =  E(D  +  S)  =  B  +  B**P  ,  (2.1) 

where B  =  E(D ) .  If  we  use  d;,  s;,  and  y;  (i=l,...,N;  N  is  the  total  number  of  pixels  in 
the  image)  to  represent  the  elements  in  the  matrices  D,  S  and  Y,  then  Poisson  noise 
model  is  as  follows: 

dj  \  B  ~  Poisson{bi ) 

Sj\B  ~  Poisson{(B  *  *P).)  ^  2) 

v,  =  di  +  S'  I  B  ~  Poisson(bl  +  (B  **P).) 


where  d;  and  s;  (i=l,...,N)  given  B  are  mutually  independent.  The  purpose  of  using 
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N 

(B**P)j  to  represent  'Y  b:pl:  is  two-folds:  1)  (B**P);  is  more  straightforward  than 

7=1 

N 

YbjPij  and  2)  it  is  a  reminder  to  us  that  two-dimensional  convolution  has  a  fast 

7=1 

implementation  in  Fourier  domain.  We  are  interested  in  estimating  B  =  {bi  ; 
i=l,...,N}. 

In  [26],  a  detailed  derivation  of  MLE  estimators  of  B  through  Expectation 
Maximization  was  provided.  For  conciseness,  only  the  final  iterative  equation  is 
shown  here: 


b 


(«+ 1) 
k 


(«) 


yk 


bf  +(BW  **P) 


(»)  *  ; 


(2.3) 


2.3  Gaussian  Noise  Model 

The  radiation  is  intrinsically  related  to  the  number  of  photons  via  a  conversion  factor, 
which  is  a  function  of  the  photon  energy.  Even  if  an  ideal  monochromatic  x-ray  beam 
is  available,  i.e.  all  the  photons  from  the  x-ray  source  carry  the  same  energy,  the 
detected  primary  photons  have  the  same  energy  but  the  detected  scatter  photons  vary 
in  their  energy  levels,  thus  causing  different  radiation  or  exposure  conversion 
efficiency.  This  issue  will  be  even  more  apparent  when  in  practice  a  polychromatic 
x-ray  beam  as  what’s  shown  in  Figure  1.1  is  usually  used.  In  this  case  even  the 
primary  photons  will  take  on  different  energy  levels.  The  radiation  or  exposure  will 
be  related  not  only  to  the  number  of  photons,  but  also  to  their  individual  energies  and 
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the  energy-dependent  conversion  factors. 

As  will  be  discussed  in  the  next  chapter,  our  experimental  data  show  that  primary 
d;  and  scatter  radiation  Sj  can  not  be  directly  modeled  as  Poisson  distribution.  By 
contrast,  the  data  approximately  follow  Gaussian  distribution.  Thus,  a  Gaussian  noise 
model  is  proposed  as  follows: 

dj  I  B,<Jn  ~  Gaussian(bj,(7n~ ) 

s:  I  B,cri2 2  ~  Gaussian((B**P)i,cri22)  (2.4) 

yf  =di  +Sj  I  B,<jn2,<Ji22  ~  Gaussian(bi  +  (B**P)i,crn 2  +crj22) 

where  di,  s;,  yi  and  b;  have  the  same  meaning  as  those  in  Poisson  noise  model  (in 

2  2 

block  (2.2)).  In  addition,  a  n  anda  &  represent  the  variance  of  the  primary 


radiation  and  the  variance  of  the  scatter  radiation  in  each  pixel  i. 

Due  to  the  convolution  operation,  the  estimation  of  B  =  {bj;  i=l,...,N}  directly 
from  Y  does  not  have  a  simple  analytic  form.  The  MLE  of  B  is  thus  derived  through 
the  EM  algorithm  as  follows. 

Treat  the  measured  Y  =  {y;,  i=l,...,N}  as  an  incomplete  dataset,  and  unobserved 

(D,S)  =  { (di,Si),  i=l, _ ,N}  as  a  complete  dataset.  The  d;  ’s  and  S;  ’s  given  B  are 

mutually  independent,  therefore  the  complete  data  likelihood  is: 


pc{D,S\BAcr2,C7i2-i=l...m>fl-T^ 


-idj-bjfllcTji2 


j= 1  J2 7ZG-, 


(2.5) 

Assuming  {aii  ,  a;2  ;  i=l,...,N}  are  known,  we  will  get  the  complete  data  log 


likelihood  by  taking  the  logarithm  on  both  sides, 
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L.  (B I  D,S)  =  i\-(d,  - bjf /  2a/  -  ( Sj-(b *  *p)/  /  2a/  -  log  /xa/  -  log  /xa/]  ■ 

M 

(2.6) 

The  EM  algorithm  is  comprised  of  two  steps:  one  is  the  E-step  where  the 
expectation  of  the  complete  data  log  likelihood  with  respect  to  the  present  estimate 
of  B  is  computed;  and  the  other  is  the  M-step  where  a  new  estimate  of  B  is  obtained 
which  will  maximize  the  computed  expectation  in  the  E-step. 

Firstly,  let  us  consider  the  E-step: 

Q(B\B(n) ) 

=  E[LC(B\D,S)\Y,BW] 

N 

=  £ {-(b/  -  2d /)bJ ) / 2c>  / 2  -  [{B  *  *P) .  )2  -  2 Sj{n) ( B**P ) . } / 2a/ 

i= i 

+  terms  ■  independent  -  of  ■  B), 

(2.7) 

d/  =E[di\Y,B(n)] 

where  .  (2.8) 

s/  =E[Sj\Y,Bln) 1 

Secondly,  consider  the  M-step  to  find  B(n+1)  that  will  maximize  Q(BIB<n)): 

SQ(B\  Bin>)  _  o  _  _(2bi  _  2dkm)l2aP  -  ^[2(fl  »  *P);  pf  -2s"' ,]l2vp . 
Solving  the  above  equation  for  bk  gives 

b/+l)  =  d/n) -cj/jZPp  -[(5("+1)  **p)j -s/Vct/.  (2.9) 

i= i 

Using  B<n)  to  approximate  B<n+1)  in  the  right  hand  side,  we  get: 
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(2.10) 


b™  =dkw-akl2^pjk  •[( B **P)j -Sjin)l/crj22. 
j= i 

As  a  good  estimate  of  the  primary  image  is  formed,  ( B[n)  *  *P) .  -  s;.(n)  «  0 ,  then, 


b[n+V)  =d[n).  (2.11) 

The  same  apparent  form  was  obtained  for  Poisson  noise  model  in  (261.  But  due  to  the 
different  statistical  models,  the  actual  forms  of  dklnl  are  different  and  so  do  the 
iterative  formula  of  bk.  We  will  get  the  iterative  formula  of  bk  for  the  Gaussian  noise 
model  through  the  following  theorem. 

Theorem  2.1  Let  X  ~  Gaussian(px,ox2),  Y  ~  Gaussian)  py,oy2),  independent.  Let 
Z=X+Y  be  the  third  random  variable.  We  know  that  Z  follows  a  Gaussian 
distribution  with  mean  of  pz  =(px+py)  and  variance  of  oz"  =(gx  +ay  ).  It  can  be 

proved  that  the  conditional  distribution  of  XIX+Y  i.e.  XIZ  is 

2  2  2  2  2 

<7  <T  U  cr  <7 

Gaussian((-^j  Z  +  / ux - ^  juy ),  — — —) . 

Gz  Gz 

Proof:  According  to  the  definition  of  conditional  distribution: 


=  p(X,Z)  =  p(Z\X)p(X) 
P(Z )  i HZ) 


(2.12) 


Note  that 

p(Z\X)  =  p(X  +Y\X)  =  p(Y)  =  p(Z-X).  (2.13) 

Therefore,  equation  (2.12)  becomes: 

P(Z-X)P(X) 

P(Z) 
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[x-(%Z+%iux-%iu,)f 


2  2 
<TX  CTV 

2(— — 2^“ ) 
a," 


2  2 

(J  <J 

2  K  x  / 

'  cr 


~  Gaussian((-^Z  H — V  // 


2  2 

cr  cr 

X  J  ' 


2  2  ^  x  2  r~y* 

G  _  G_  G_ 


Because  the  primary  and  scatter  radiation  of  each  pixel  given  B  is  independent  of 
those  of  other  pixels,  dk(n)  =  E[dk \Y,B{n}]  =  E[dk  \  yk  =  dk  +  sk,B{,,)]  .  Consider 
random  variables  in  Theorem  2.1  to  be  X=dk,  Y=Sk  and  Z=yk,  therefore, 


dkn)  -<r n(Bin)  +o-,22). 


(2.14) 


Then  equation  (2.11)  combines  with  equation  (2.14)  to  give  the  following  updating 
equation: 


b'r"  =i>,“  +W„  [v.  -<V  +  (B<"> 


(2.15) 


where 


Wk  +^22)- 


(2.16) 
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2.4  Incorporation  of  Gibbs  Prior 


MLEM  is  known  to  have  adverse  effect  on  high  frequency  image  noise.  To 
overcome  this,  some  constraints  can  be  put  on  the  noise  level  within  the  estimated  B, 
or  in  other  words,  the  prior  information  about  B  is  provided.  By  Bayes’s  Rule, 
P(B\Y)kP(Y\B)P(B),  (2.17) 

where  p(B)  is  the  prior  joint  distribution  of  B={b;;  i=l,...,N},  p(YIB)  is  equal  to  the 
likelihood  of  B,  and  p(BIY)  is  the  posterior  joint  distribution  of  B  given  measured 
pixel  values  Y={yg  i=l,...,N}. 

We  assume  B  is  a  Markov  random  process;  it  therefore  follows  a  Gibbs 
distribution: 


K 


(2.18) 


where  K  is  a  normalizing  factor  which  is  independent  of  B,  U(B)  is  the  energy 
function,  and  P  is  a  free  parameter  adjusting  the  relative  weight  of  this  prior  on  the 
maximum  a  posteriori  (MAP)  estimator  of  B.  When  P  is  approaching  infinity,  the 
MAP  of  B  approaches  MLE  of  B. 

The  energy  function  is  the  sum  of  the  potential  function,  i.e., 

U{B)  =  YVC{B ),  (2.19) 

ceC 

where  C  is  the  set  comprised  of  all  cliques  in  the  image.  One  clique  is  defined  as  a 
set  of  pixels  where  each  one  is  a  neighbor  of  all  the  others  in  the  same  clique.  In  the 
thesis,  the  Gibbs  prior  is  defined  over  a  2nd  -order  neighborhood  system  (for  each 
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pixel  its  north,  south,  east,  west  neighboring  pixels  plus  its  four  diagonal  neighboring 
pixels)  with  each  clique  comprising  of  two  neighboring  pixels.  There  are  many  forms 
of  the  potential  function  Vc(b).  The  one  we  pick  up  is  adaptive  to  discontinuity 
[35-37]: 


vabr-bj}) 


(b.-bj)2 

sS  +  tb'-bj)2’ 


(2.20) 


where  i  and  j  are  the  neighboring  pixels  within  the  clique  i~j.  The  bi  and  bj  represent 
their  intensities.  5C  is  an  adjustable  parameter  to  regulate  the  cut-off  frequency  of  the 
noise  in  the  image. 


2.5  Image  Acquisition 

A  Siemens  prototype  digital  mammography  system  (Mammomat  Novation  )  with 
70  pm  isotropic  resolution  was  used  for  image  acquisition.  Uniform  breast  phantoms 
(CIRS,  Inc.,  Norfolk,  VA)  were  imaged  with  the  x-ray  beam  generated  by  28kVp  and 
Mo/Mo  target/filter  combination.  The  phantoms  are  radiographically  equivalent  to  a 
compressed  breast  of  4cm  in  thickness  and  50%  in  glandular  tissue  density.  At  the 
center  of  the  phantom  there  is  a  square  dent,  which  mimics  a  high-contrast  lesion  in 
the  digital  mammography  images.  The  images  used  in  this  thesis  were  acquired 
without  an  anti-scatter  grid.  Some  of  them  have  a  beam  stop  (i.e.,  lead  discs  with 
3mm  diameter)  array  superimposed  on  the  breast  phantoms.  The  beam  stop  method 
is  a  standard  technique  to  measure  the  scatter  radiation.  Because  lead  discs  absorb  all 
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the  primary  radiation,  only  scatter  radiation  can  arrive  behind  them.  Figure  2.2  shows 


one  such  image. 


Figure  2.2:  A  sample  breast  phantom  image  taken  with  the  beam  stop  array 
superimposed.  The  bright  square  mimics  a  high  density  lesion,  based  on  which  the 
contrast  and  CNR  values  are  obtained. 


The  image  can  then  be  fed  into  the  algorithms  for  processing.  The  effect  of 
processing  is  evaluated  through  various  metrics,  which  will  be  discussed  in  the 
following  subsection. 


2.6  Image  Analysis  Metrics 


The  primary  purpose  of  the  algorithms  is  to  estimate  the  expectation  of  primary 
radiation.  Its  effect  is  measured  by  the  residual  scatter  fraction  (RSF).  At  the  same 
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time,  it  is  desirable  that  the  contrast-to-noise  ratio  (CNR)  will  be  constrained  or  even 


improved  after  image  processing.  In  addition,  the  effect  of  the  algorithms  on  spatial 
resolution  of  the  image  has  to  be  carefully  monitored.  In  the  following,  we  will  give 
the  definition  of  each  of  these  metrics  and  how  they  are  measured  in  this  thesis. 


2.6.1  Residual  Scatter  Fraction 


Scatter  fraction  (SF)  is  defined  as  the  ratio  of  the  scatter  radiation  to  the  total 
radiation.  Residual  scatter  fraction  (RSF)  is  a  quantity  used  to  indicate  how  much  of 
the  scatter  radiation  remains  after  applying  the  scatter  compensation  algorithm. 

For  the  given  imaging  technique,  two  sets  of  images  of  the  phantom  were 
obtained.  One  is  taken  without  a  beam  stop  array,  and  the  other  is  taken  with  the 
beam  stop  array.  The  signals  behind  beam  stops  (lead  discs)  are  the  scatter  radiation. 
The  total  radiation,  which  is  the  sum  of  primary  radiation  and  the  scatter  radiation, 
will  reach  the  region  without  the  beam  stops.  Thus  the  measured  primary  radiation 
(Pmeasured)  is  calculated  by  subtracting  the  mean  radiation  of  a  region-of-interest  (ROI) 
behind  a  beam  stop  from  the  mean  of  the  same  ROI  location  without  a  beam  stop.  In 
the  image  processed  for  scatter  compensation,  the  mean  of  total  radiation  (T)  in  the 
same  ROI  location  (Testimated)  is  the  sum  of  the  residual  scatter  radiation  and  the 
primary  radiation.  Then 


RSF 


T  -  P 

estimated  measured 

T 

estimated 


(2.21) 
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2.6.2  Contrast,  Noise  and  CNR 


The  contrast  is  defined  as  the  ratio  of  the  difference  between  the  mean  value  of  the 
lesion  (Tiesion)  and  that  of  the  background  (Tbackground)  to  the  mean  of  the  background. 
That  is, 

T  -T 

Contrast  = - - - .  (2.22) 

T 

background 

The  noise  is  derived  by  dividing  the  standard  deviation  (STDbackground)  to  the 
mean  (Tbackground)  of  the  background: 

STD 

Noise  = - ba*g™nd  .  (2.23) 

T 

background 


Contrast-to-noise  ratio  is  the  ratio  of  the  contrast  to  the  noise,  i.e., 


CNR  = 


Contrast  3 'lesion 


lesion  1  background 


Noise 


STD 


(2.24) 


background 


2.6.3  Resolution 

Due  to  the  nonlinearity  of  the  algorithm,  the  metric  like  modulation  transfer  function 
(MTF)  which  is  designed  for  a  linear  system  can  not  be  used  here.  Instead,  a  test  bar 
comprised  of  alternating  bright  and  dark  lines  with  size  corresponding  to  Nyquist 
frequency  are  embedded  in  the  phantom  image. 

The  contrast  improvement  factor  (CIF),  defined  as  the  ratio  of  the  contrast  after 
image  processing  to  the  initial  contrast,  is  obtained  for  the  test  bar  with  various  initial 
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contrast.  If  CIF  is  not  less  than  1,  no  resolution  is  lost.  Otherwise,  resolution  is 


degraded.  The  minimal  initial  contrast  that  the  test  bar  can  take  on  with  CIF  no  less 
than  1  is  recorded  as  an  indication  of  the  effect  of  the  image  processing  on 
resolution. 
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Chapter  3 

Results  and  Discussion 


3.1  Normality  Check  of  the  Data 

In  our  new  model,  primary,  scatter  and  total  radiations  of  each  pixel  in  the  projection 
mammography  images  are  assumed  to  follow  Gaussian  distributions.  To  check 
whether  Gaussian  distribution  is  a  good  approximation  to  the  real  data,  we  analyze  a 
uniform  ROI  outside  of  a  beam  stop  (ROI1)  and  a  uniform  ROI  behind  a  beam  stop 
(ROI2)  in  an  image  acquired  without  an  anti-scatter  grid. 

ROI1  is  a  square  region  with  201x201  pixels,  and  its  histogram  is  plotted  in 
Figure  3.1.  Visually,  it  follows  approximately  a  Gaussian  distribution.  For  further 
evaluation,  the  empirical  quantile-quantile  plot  of  the  data  with  respect  to  a  standard 
Gaussian  distribution  is  drawn  in  Figure  3.2.  The  quantile  of  the  data  has  a  nice 
linear  relationship  with  the  quantile  of  the  standard  Gaussian  distribution,  indicating 
that  the  total  radiation  can  be  well  represented  by  a  Gaussian  distribution. 

As  stated  before,  the  exposure  of  the  area  behind  a  beam  stop  is  due  to  the 
scattered  radiation  from  the  neighboring  regions.  Figure  3.3  shows  the  profile  of 
radiation  along  a  line  through  the  center  of  a  beam  stop.  It  has  a  nice  flat  profile  for 
the  scatter  radiation.  The  circular  ROI2  is  selected  which  has  a  total  number  of  441 
pixels.  Its  histogram  and  quantile-quantile  plot  with  respect  to  the  standard  Gaussian 
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distribution  are  shown  in  Figure  3.4  and  Figure  3.5  respectively.  Figure  3.5  shows 
that  scatter  radiation  is  also  approximately  Gaussian  distribution. 

It  is  obvious  that  both  total  radiation  shown  in  Figure  3.1  and  scatter  radiation 
shown  in  Figure  3.4  can  not  be  directly  modeled  as  Poisson  distributions,  since  1)  the 
radiation  does  not  take  on  discrete  integer  values  only;  and  2)  the  variance  is  much 
smaller  than  the  mean,  whereas  a  Poisson  distribution  has  an  equal  variance  and 
mean.  Therefore,  from  the  modeling  perspective,  the  Poisson  noise  model  (shown  in 
block  (2.2))  is  problematic  especially  when  a  computation  method  other  than  EM 
algorithm  is  used.  Luckily,  if  the  model  is  modified  by  adding  a  latent  data,  then  the 
iterative  equation  (equation  (2.3))  derived  from  EM  algorithm  can  be  approximately 
true  under  certain  assumptions. 


Figure  3.1:  The  histogram  of  data  from  a  uniform  region-of-interest  with  total  of 
201x201  pixels  in  an  image  acquired  without  an  anti-scatter  grid.  The  data  is  seen  to 
be  approximately  Gaussian  distribution. 
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Figure  3.2:  The  quantile-quantile  plot  of  the  same  data  as  in  Figure  3.1  with  respect 
to  the  standard  Gaussian  distribution  with  mean  of  0  and  standard  deviation  of  1 .  The 
data  fits  well  with  the  Gaussian  distribution. 


Figure  3.3:  One-dimensional  profile  through  the  center  of  a  beam  stop  (or  lead 
disc). 
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Figure  3.4:  The  histogram  of  data  from  a  circular  region-of-interest  with  total  of  441 
pixels  behind  a  beam  stop  in  an  image  acquired  without  an  anti-scatter  grid.  The  data 
is  seen  to  be  approximately  Gaussian  distribution. 
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Figure  3.5:  The  quantile-quantile  plot  of  the  same  data  as  in  Figure3.4  with  respect  to 
the  standard  Gaussian  distribution  with  mean  of  0  and  standard  deviation  of  1. 


27 


Except  an  outlier  (data  value  is  about  1.7mR),  the  rest  data  fits  well  with  the 
Gaussian  distribution. 


3.2  Latent  Data  Augmentation  Scheme  for  Poisson  Model 


Assume  the  exposure  or  radiation  is  intrinsically  proportional  to  the  number  of 
photons  that  produce  the  radiation  through  a  constant  C  (E),  which  is  dependent  on 
the  energies  of  photons  E.  For  di,  Si  and  yh  the  corresponding  number  of  photons  are 
nd;,  ns;,  and  nyj.  The  expected  number  of  photons  for  ndj  is 
nb;.  The  following  model  is  valid: 


ndt  I  NB  ~  Poisson(nbj) 

nSj  I  NB  ~  Poisson((B**P)j ) 

nVj  =  ndj  +  ns j  I  NB  ~  Poissoninbj  +  (NB**P)j) 


(3.1) 


The  updating  equation  for  nb;  is  the  same  as  equation  (2.3): 


nbk"+l  1  =  nbk 


(«) 


nyk 


nbi"’  +(NBin)  **P). 


(3.2) 


If  by  simplification,  assume  a  single  energy  independent  conversion  constant  C 
instead  of  a  set  of  C(E)  exists,  then 


bk  =  E(dk )  =  C  •  E(ndk )  =  C  ■  nbk 
yk  =C-nyk 


(3.3) 


And  then, 


,(«+!)  _  /, 
'b  , 


yk 


bk  +  (Bw  **P). 


(3.4) 


which  is  the  same  updating  equation  for  bk  as  equation  (2.3).  In  other  words,  as  long 
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as  a  single  conversion  constant  rather  than  a  set  of  energy  dependent  conversion 
constants  C(E)  exists  between  the  radiation  and  the  number  of  photons,  The  equation 
(2.3)  or  (3.3)  derived  from  EM  algorithm  is  scale-invariant  regardless  of  the  value  of 
latent  data  C. 

3.3  Comparison  between  the  Poisson  and  Gaussian  Models 

Gaussian  noise  model  is  sounder  than  Poisson  noise  model  for  the  direct  modeling  of 
radiations.  Poisson  model  won’t  give  an  accurate  answer  if  computation  methods 
which  rely  heavily  on  the  accuracy  of  the  model  such  as  Gibbs  sampling  are  adopted. 
However,  when  EM  algorithm  is  used,  the  iterating  equation  obtained  based  on 
Poisson  assumption  can  roughly  be  used  to  update  the  expectation  of  primary 
radiation.  This  is  based  on  a  simplification  that  a  single  conversion  factor  is  valid  for 
all  radiation  levels,  which  is  not  true  in  reality.  By  contrast,  Gaussian  noise  model 
allows  for  the  different  conversion  factors  for  different  energies. 

For  the  uniform  breast  phantom,  empirical  value  of  Wi  (i=l,...,N)  in  the 
Gaussian  model  as  shown  in  equation  (2.16)  is  0.45.  The  empirical  optimized  scatter 
kernel  P  has  FWHM=80pixels  (i.e.,  5.6mm)  and  M=0.52.  The  images  are  processed 
to  obtain  the  MLE  estimates  from  both  models.  Also  the  MAP  estimates  are  obtained 
from  both  models  using  the  same  Gibbs  prior  with  delta  of  0.10. 

As  a  convention,  in  all  the  figures  shown  hereafter,  the  values  for  iteration  0  are 
the  values  measured  on  the  original  image  without  any  processing. 
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Figure  3.6  shows  the  plots  of  RSF  of  MLE  and  MAP  estimates  from  both  models 
as  a  function  of  iteration  number.  All  estimators  successfully  reduce  the  scatter 
radiation  in  the  processed  image  such  that  RSF  drops  with  iteration.  In  addition,  they 
all  asymptotically  converge  to  the  same  value  of  0.019  from  the  original  scatter 
fraction  of  0.354.  Note  that  estimators  based  on  Poisson  noise  model  have  a  slightly 
faster  convergence  rate  than  those  based  on  Gaussian  noise  model.  For  example,  at 
iteration  2,  the  RSF  of  estimators  from  Poisson  model  already  drops  to  0.019, 
whereas  that  from  Gaussian  model  is  0.047. 


— * — MLE  -  Poisson  Model  MLE  -  Gaussian  Model 

— *—  MAP  -  Poisson  Model  MAP  -  Gaussian  Model 


Figure  3.6:  RSF  vs.  iteration  plots  for  MLE  and  MAP  estimator  of  b  from  both  the 
Poisson  noise  model  and  the  Gaussian  noise  model.  The  magnitude  of  scatter  kernel 
is  0.52,  which  is  same  as  the  measured  scatter- to-primary  ratio  (SPR).  Therefore, 
RSF  all  drops  close  to  zero,  meaning  almost  complete  scatter  compensation. 


Figure  3.7  to  Figure  3.9  illustrate  how  noise,  contrast  and  CNR  individually 
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change  with  iteration  numbers.  Let’s  look  at  Figure  3.7  first.  At  iteration  16,  the 
MLE  estimators  of  both  Poisson  model  and  Gaussian  model  increase  the  noise  from 
the  original  0.021  to  0.031,  which  corresponds  to  a  47.6%  of  increase.  The  MAP 
estimator  from  Poisson  model  keeps  the  noise  at  a  roughly  same  level  as  the  original 
image,  whereas  the  MAP  estimator  from  Gaussian  model  decreases  the  noise  by 
5.6%. 

Figure  3.8  show  that  all  four  estimators  increase  the  contrast  similarly.  At 
iteration  16,  they  all  increase  the  contrast  by  26.8%. 

As  shown  in  Figure  3.9,  the  initial  CNR  of  the  lesion  is  47.28.  After  processing 
by  both  MLE  methods,  the  CNR  of  the  lesion  drops  to  39.13,  which  is  equivalent  to 
a  17.2%  change.  By  contrast,  the  MAP  estimators  from  Poisson  and  Gaussian 
models  successfully  increase  CNR  by  18.7%  and  34.2%  respectively. 

Table  3.1  shows  the  resolution  results  for  the  MAP  estimates  from  the  Poisson 
noise  model  and  the  Gaussian  noise  model.  Both  models  can  retain  resolution  for 
initial  contrast  greater  than  2%.  Poisson  model  retains  the  resolution  slightly  better 
than  Gaussian  model.  It  is  not  shown  in  Table  3.1  that  MLE  estimates  from  both 
models  retain  the  resolution  at  all  initial  contrast  levels. 

In  summary,  both  MLE  and  MAP  estimators  work  equally  well  in  reducing  the 
scatter  radiation,  while  MAP  estimators  works  better  than  their  MLE  counterparts  in 
improving  or  constraining  CNR  without  general  loss  of  resolution.  MAP  estimator 
based  on  Gaussian  model  has  a  better  performance  in  improving  or  constraining 
CNR  than  the  MAP  estimator  based  on  Poisson  model. 
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Figure  3.7:  The  noise  vs.  iteration  plots  for  MLE  and  MAP  estimator  of  b  from  both 
the  Poisson  noise  model  and  the  Gaussian  noise  model. 
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Figure  3.8:  The  contrast  noise  vs.  iteration  plots  for  MLE  and  MAP  estimator  of  b 
from  both  the  Poisson  noise  model  and  the  Gaussian  noise  model. 


Figure  3.9:  The  CNR  vs.  iteration  plots  for  MLE  and  MAP  estimator  of  b  from  both 
the  Poisson  noise  model  and  the  Gaussian  noise  model. 


Poisson  Model 

Gaussian  Model 

Minimum  initial  contrast  that 

is  retainable  during  processing 

1.8% 

2.0% 

Table  3.1:  The  resolution  results  for  Poisson  model  and  Gaussian  model.  The  square 
wave  function  with  Nyquist  frequency  is  used  as  the  test  object.  For  various  initial 
contrasts,  the  corresponding  contrast-improvement-factor  (CIF)  is  computed  at 
iteration  16.  CIF  no  less  than  1  is  used  as  the  criterion  for  retaining  the  spatial 
resolution.  What  is  reported  here  is  the  minimal  initial  contrast  that  has  CIF  no  less 
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than  1. 


3.4  Further  Evaluation  of  Gaussian  Noise  Model 

3.3.1  Effect  of  the  Magnitude  of  the  Scatter  Kernel 

The  Magnitude  of  the  Scatter  Kernel  M,  which  is  the  same  as  the  area  under  the 
curve,  is  used  to  model  the  scatter-to-primary  ratio  (SPR).  Using  the  specified 
technique,  the  measured  SPR  for  the  phantom  is  0.52.  As  is  shown  in  Figure  3.10, 
when  M  is  specified  as  0.52,  the  RSF  drops  rapidly  from  the  initial  value  to  a  value 
close  to  zero,  meaning  a  satisfactory  scatter  compensation  effect.  When  M  is  less 
than  0.52,  the  scatter  radiation  is  partially  compensated.  Specifically,  when  M  equals 
zero,  no  scatter  compensation  is  made.  When  M  is  larger  than  0.52,  the  scatter 
radiation  is  over-compensated,  i.e.,  RSF  is  less  than  zero. 
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Figure  3.10:  The  RSF  vs.  iteration  curves  for  the  magnitude  (M)  of  the  scatter  kernel 
ranging  from  0.0  to  0.65.  The  measured  SPR  value  using  the  beam  stop  technique  is 
0.52.  M  of  0.0  represents  no  scatter  compensation.  M  of  0.20  and  0.40  represents 
partial  scatter  compensation.  And  M  of  0.65  overcompensate  the  scatter  radiation  in 
the  image. 


Figure  3.11  -3.13  illustrate  how  different  magnitude  of  the  scatter  kernel  affects 
the  noise,  contrast  and  CNR  respectively.  Overall,  these  three  metrics  change 
monotonically  with  respect  to  the  magnitude  M.  More  specifically,  when  M  gets 
larger,  both  the  noise  and  the  contrast  become  larger,  while  the  CNR  gets  smaller. 

For  all  the  M  values  investigated,  the  CNR  at  iteration  16  is  larger  than  the 
original  value.  When  M  is  equal  to  0,  the  CNR  improves  by  as  large  as  74.7%;  even 
for  M  of  0.65,  CNR  improves  by  24.3%. 

Table  3.2  gives  the  resolution  results  for  the  Gaussian  model  with  various 
magnitude  of  scatter  kernel.  Note  that  for  a  magnitude  of  zero,  there  will  always  be 
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resolution  loss.  For  the  rest  magnitude  values,  the  resolution  performances  are 
similar. 

If  the  scatter  compensation  is  the  major  concern,  then  the  magnitude  should  be 
chosen  as  close  to  the  actual  SPR  value  as  possible.  The  image  processed  in  this 
setting  will  improve  or  constrain  CNR  without  general  loss  of  resolution.  If  more 
noise  reduction  or  more  CNR  improvement  is  desired,  then  a  smaller  magnitude  can 
be  selected,  at  the  expense  of  partial  scatter  compensation.  The  magnitude  of  zero, 
however,  is  not  a  good  choice  because  it  reduces  the  noise  and  increases  the  CNR  at 
the  expense  of  all  spatial  resolution  as  well  as  no  scatter  compensation. 
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Figure  3.11:  The  noise  vs.  iteration  curves  for  M  ranging  from  0.0  to  0.65.  At  each  M 
level,  the  percentage  noise  reduces  asymptotically.  The  smaller  the  M  value  is,  the 
more  is  the  percentage  noise  reduced  from  the  initial  value  of  0.021. 
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Figure  3.12:  The  contrast  vs.  iteration  curves  for  M  ranging  from  0.0  to  0.65.  At  each 
magnitude  level,  the  contrast  increases  asymptotically.  The  larger  the  magnitude  is, 
the  more  is  the  contrast  increased  from  the  initial  value  of  0.97. 
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Figure  3.13:  The  CNR  vs.  iteration  curves  for  M  ranging  from  0.0  to  0.65.  At  each 
magnitude  level,  the  contrast  increases  asymptotically.  The  smaller  the  magnitude  is, 
the  more  is  the  contrast  increased. 


Magnitude 

0.0 

0.2 

0.4 

0.52 

0.65 

Minimal 

initial  contrast 

2.7 % 

2.2% 

2.0% 

2.0% 

Table  3.2:  The  resolution  results  for  the  Gaussian  noise  model  with  different 
magnitude  of  scatter  kernel. 

3.3.2  Effect  of  the  Delta  in  the  Gibbs  Prior 

As  mentioned  in  the  subsection  2.4,  the  delta  (5)  in  the  potential  function  of  the 
Gibbs  prior  can  be  considered  as  a  factor  controlling  cut  off  frequency  in  the 
processed  image.  The  results  shown  in  the  previous  sections  are  for  8=0.10.  Now 
let’s  consider  5  of  0.2  and  0.05  to  see  how  plots  for  the  image  quality  metrics  change. 
As  expected,  different  5  values  do  not  change  the  RSF  and  the  contrast  plots,  so  their 
plots  are  not  shown.  Figure  3.14,  Figure  3.15  and  Table  3.3  give  the  noise,  CNR 
plots  and  the  resolution  result  for  delta  of  0.2.  Figure  3.16,  Figure  3.17  and  Table  3.4 
show  the  corresponding  results  for  delta  of  0.05.  By  comparing  these  results  to 
Figure  3.7,  Figure  3.9  and  Table  3.1  for  delta  of  0.1,  the  trend  emerges:  the  larger  5 
is,  the  more  the  image  will  be  smoothed  at  the  expense  of  slightly  more  resolution 
loss.  It  is  understandable  since  more  smoothed  the  final  image  is  (i.e.,  less  noise  in 
the  final  image),  more  likely  the  details  will  lose.  But  still  the  final  image  retains  a 
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reasonably  good  spatial  resolution  (for  test  bar  of  smallest  possible  size 
(corresponding  to  Nyquist  frequency),  all  initial  contrasts  of  3%  or  larger  is  retained). 
When  delta  is  0.2,  CNR  in  the  MAP  estimator  based  on  the  Gaussian  model 
improves  by  as  large  as  120%  without  general  loss  of  resolution. 


Figure  3.14:  The  noise  vs.  iteration  plots  for  MLE  and  MAP  estimators  based  on 
delta  of  0.2.  The  MAP  estimators  decrease  the  noise  more  than  their  counterparts 
based  on  delta=0.1  as  shown  in  Figure  3.6. 
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Figure  3.15:  The  CNR  vs.  iteration  plots  for  MLE  and  MAP  estimators  based  on 
delta  of  0.2.  The  MAP  estimators  increase  the  CNR  more  than  their  counterparts 
based  on  delta=0.1  as  shown  in  Figure  3.8.  Also,  the  MAP  estimator  from  Gaussian 
model  improves  CNR  more  than  the  one  from  Poisson  model. 


Poisson  Model 

Gaussian  Model 

Minimum  initial  contrast  that 

is  retainable  during  processing 

3.0% 

3.1% 

Table  3.3:  The  resolution  results  for  Poisson  model  and  Gaussian  model  when 
delta=0.2. 
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Figure  3.16:  The  noise  vs.  iteration  plots  for  MLE  and  MAP  estimators  based  on 
delta  of  0.05.  The  MAP  estimators  decrease  the  noise  less  than  their  counterparts 
based  on  delta=0.1  as  shown  in  Figure  3.6. 
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Figure  3.17:  The  CNR  vs.  iteration  plots  for  MLE  and  MAP  estimators  based  on 
delta  of  0.05.  In  this  case,  the  MAP  estimators  performs  better  than  the  MLE 
estimators  in  constraining  CNR,  but  performs  worse  than  their  counterparts  based  on 
delta=0.1  as  shown  in  Figure  3.8.  But  the  MAP  estimator  from  Gaussian  model  is 
still  slightly  better  than  the  one  from  Poisson  model. 


Poisson  Model 

Gaussian  Model 

Minimum  initial  contrast  that 

is  retainable  during  processing 

1.2% 

1.6% 

Table  3.4:  The  resolution  results  for  Poisson  model  and  Gaussian  model  when 
delta=0.05. 
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Chapter  4 
Conclusion 


By  checking  the  experimental  data,  it  was  found  that  the  Poisson  noise  model  for 
scatter  compensation  in  the  literature  [26]  can  not  account  for  the  radiations 
(including  the  primary,  scatter  and  total  radiations)  directly.  It  will  lead  to  an 
erroneous  result  for  the  estimation  of  the  expected  values  of  primary  radiation  if  a 
computation  method  like  Gibbs  sampling  is  used.  Luckily,  due  to  the 
scaling-invariant  property  of  EM  algorithm  with  an  approximation  that  a  single 
factor  rather  than  a  set  of  energy  dependent  ones  exists  between  the  conversion  of 
radiation  and  the  corresponding  number  of  photon,  the  updating  equation  derived 
from  the  old  model  can  still  be  useful. 

The  histograms  of  radiation  data  indicate  that  they  might  be  modeled  by  a 
different  distribution  like  Gaussian.  The  quantile-quantile  plots  of  the  data  with 
respect  to  the  standard  Gaussian  distribution  show  that  Gaussian  noise  model  can  be 
reasonably  assumed.  The  EM  algorithm  based  on  this  new  model  is  derived  and 
implemented.  A  MAP  algorithm  by  incorporating  a  Gibbs  prior  is  also  implemented 
for  better  CNR  in  the  processed  images. 

The  MLE  and  MAP  estimators  from  the  Gaussian  noise  model  are  compared 
with  their  counterparts  based  on  the  Poisson  noise  model.  Results  show  that  MAP 
estimators  from  both  models  have  better  CNR  performance  than  MLE  ones  without  a 
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significant  loss  of  resolution.  In  addition,  the  MAP  estimator  from  Gaussian  model 
performs  better  than  the  one  from  Poisson  model  in  CNR  improvement. 

Further  evaluation  of  MAP  estimators  from  Gaussian  model  shows  that  both  the 
magnitude  of  the  scatter  kernel  and  the  delta  in  the  Gibbs  prior  can  be  used  to  adjust 
the  noise  and  CNR  level  in  the  processed  images.  The  delta  in  the  Gibbs  prior  acts  as 
a  major  tuner  of  CNR  without  affecting  RSF,  whereas  the  magnitude  of  the  scatter 
kernel  acts  as  a  fine  tuner  of  CNR.  Changing  the  magnitude  of  the  scatter  kernel  will 
also  affect  the  scatter  compensation  level.  There  is  a  general  tradeoff  between  the 
CNR  improvement  and  resolution  reservation.  Fortunately,  for  the  largest  CNR 
improvement  (2.2  times  the  original  CNR),  the  resolution  is  still  reasonably  well 
reserved. 
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